#!/usr/bin/python
# GM_Locations.py

import math
import numpy

VariableNames = ['dlocs', 'delds', 'PtScaled', 'llocs', 'dells', 'PtScalela', 'blocs', 'delbs', 'PtScaleba']

def Locations(simpars, phypars):
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd

    # Earth-centric coords, 1d
    # d
    dlocs = numpy.arange(dres) + 0.5
    delds = numpy.ones(dres, dtype=float)
    # l
    llocs = numpy.arange(lres) + 0.5
    dells = numpy.ones(lres, dtype=float)
    # b
    blocs = numpy.arange(bres) + 0.5
    delbs = numpy.ones(bres, dtype=float)

    # 1d
    # [d] d is distance from earth in pc
    d = dlocs*PtScaled
    deltad = delds*PtScaled
#    dd = numpy.array(deltad)
    # [l] la is galactic longitude in radians instead of degrees.  l is zero when pointed towards the galactic center, pi when pointed away.
    la = llocs*PtScalela
    deltala = dells*PtScalela
#    dla = numpy.array(deltala)
    # [b] ba is galactic lattitude in radians instead of degrees.  b is zero when pointed towards the galactic plane, pi/2 when pointed up out of it, and -pi/2 when pointed down into it.
    ba = blocs*PtScaleba
    deltaba = delbs*PtScaleba
#    dba = numpy.array(deltaba)

    # Special exception for the pencilbeam case.  Must be done before higher-dimensional quantities are calculated.
    # The beam is one square degree wide - not that it matters; that width gets taken out at the end.
    if pencilBeam == True:
        # [l], radians
        lvalue = pencilBeamPointing[0]
        if (lvalue < 0) or (lvalue > 360):
            print "The value of galactic latitude l given in pencilBeamPointing must be between 0 and 360 in units of degrees, i.e. in standard notation for galactic coordinates.  Fix it.  Because this program will give you an unreliable result if you do not."
        if lvalue > 180:
            lvalue = 180 - (lvalue - 180)
        lvalue = lvalue*math.pi/180. # degrees -> radians
        la = numpy.array([lvalue])
        llocs = la/PtScalela
        deltala = numpy.array([1])*math.pi/180. # degrees -> radians
        dells = deltala/PtScalela
        lres = 1
        # [b], radians
        bvalue = pencilBeamPointing[1]
        if (bvalue < -90) or (bvalue > 90):
            print "The value of galactic longitude b given in pencilBeamPointing must be between -90 and 90 in units of degrees, i.e. in standard notation for galactic coordinates.  Fix it.  Because this program will give you an unreliable result if you do not."
        if bvalue < 0:
            bvalue = -bvalue
        bvalue = bvalue*math.pi/180. # degrees -> radians
        ba = numpy.array([bvalue])
        blocs = ba/PtScaleba
        deltaba = numpy.array([1])*math.pi/180. # degrees -> radians
        delbs = deltaba/PtScaleba
        bres = 1

    Nptsdlb = dres*lres*bres
 
    # 2d
    # [d, b]
    # z is distance out of the galactic plane in pc.
    z = numpy.multiply( d[:,numpy.newaxis], numpy.sin(ba)[numpy.newaxis,:] )
    z2 = z**2
    # d_ is the the projection of d onto the galactic plane.
    d_ = numpy.multiply( d[:,numpy.newaxis], numpy.cos(ba)[numpy.newaxis,:] )

    # [d, l, b]
    # To best align with the standard definition of galactic coordinates, I'm changing over to make my x axis be positive at l = 0, and y axis positive at l = pi/2.  Both have zero value at the galactic center.  R is distance in the galactic plane from the galactic center and r is absolute distance from galactic center.
    y = numpy.multiply(d_[:,numpy.newaxis,:], numpy.sin(la)[numpy.newaxis,:,numpy.newaxis] )
    x_ = numpy.multiply(d_[:,numpy.newaxis,:], numpy.cos(la)[numpy.newaxis,:,numpy.newaxis] )
    x = x_ - R0
    R = numpy.sqrt(x**2 + y**2)
    r2 = numpy.add(R**2, z2[:,numpy.newaxis,:])
    r = numpy.sqrt(r2)
    # phi is measured around galactic origin from positive x, and has the same sense and range as l.  It's a wee bit tricky because arctan is defined for positive x only, so I give general provisions here for negative values of x.  See notes 5/6/12 pg. 122.
    phi = numpy.arctan(y/x) # radians
    # True where x < 0
    xmask = numpy.greater(0, x)
    # Where x < 0, add pi (moving "false Q-1" to Q2 and "false Q1" to Q3; "true Q-1" i.e. Q4 can stay as is b/c all operations on it will have correct results)
    phi = phi + math.pi*xmask

    for i in phypars:
        cmd = "phypars['%s'] = %s" % (i,i)
        exec cmd
    for i in simpars:
        cmd = "simpars['%s'] = %s" % (i,i)
        exec cmd
    return phypars
